#!/usr/bin/env python

from pymesh import trigen
import matplotlib.pyplot as plt
import math

def plot_segments(m):
    for i in range(0, len(m.segments), 2):
        n1 = m.segments[i  ]-1
        n2 = m.segments[i+1]-1
        x = [m.points[2*n1],   m.points[2*n2]]
        y = [m.points[2*n1+1], m.points[2*n2+1]]
        plt.plot(x, y, "k-o")
    plt.show()

def plot_triangles(m):
    for i in range(m.numberoftriangles):
        c = get_triangle_coords(m ,i)
        x = [i[0] for i in c[0:3]]
        y = [i[1] for i in c[0:3]]
        x.append(x[0])
        y.append(y[0])
        plt.plot(x, y, "k-o")
    plt.show()

def get_triangle_nodes(m, id):
    return [m.triangles[3*id+i]-1 for i in range(m.numberofcorners)]

def get_triangle_coords(m, id):
    return [m.points[2*node:2*(node+1)] for node in get_triangle_nodes(m, id)]

def area(m):
    A = 0.
    for i in range(m.numberoftriangles):
        c = get_triangle_coords(m ,i)
        x = [j[0] for j in c[0:3]]
        y = [j[1] for j in c[0:3]]
        A += 0.5*math.fabs((x[0]-x[2])*(y[1]-y[0])-(x[0]-x[1])*(y[2]-y[0]))
    return A 

def Ix(m):
    I = 0.
    for i in range(m.numberoftriangles):
        c = get_triangle_coords(m ,i)
        x = [j[0] for j in c[0:3]]
        y = [j[1] for j in c[0:3]]
        
        dI = 0.
        for j in range(len(x)):
	        jj = (j+1)%len(x)
	        dI += (y[j]*y[j]+y[j]*y[jj]+y[jj]*y[jj])*(x[j]*y[jj]-x[jj]*y[j])
        I += (1./12.)*math.fabs(dI)
    return I 


a = trigen.MeshIO()

h = 160.
b = 160.
s = 8.
t = 13.
r = 15.
n = 150

points = [0.5*s, 0.,
          0.5*s, 0.5*h-t-r]

dtheta = 0.5*math.pi/n
for i in range(1,n+1):
    points.append(0.5*s     + r*(1.-math.cos(i*dtheta)))
    points.append(0.5*h-t-r + r*math.sin(i*dtheta))
points += [0.5*b, 0.5*h-t]
points += [0.5*b, 0.5*h]
points += [0.,    0.5*h]

temp = points[0:-2:]
for i in range(len(temp), 0, -2):
    points.append(-temp[i-2])
    points.append(+temp[i-1]) 

temp = points[2:-2:]
for i in range(len(temp), 0, -2):
    points.append(+temp[i-2])
    points.append(-temp[i-1]) 


segments = []
segments.append(1)
for i in range(len(points)/2-1):
    segments.append(i+2)
    segments.append(i+2)
segments.append(1)

# for i in range(0,len(segments),2):
#    print segments[i], segments[i+1]

a.points = points
a.segments = segments
# plot_segments(a)
# print a.points
b =trigen.MeshIO()

trigen.mesh("pqV", a, b)
#for i in range(0,len(b.points),2):
#    print i/2, b.points[i],b.points[i+1]
#for i in range(0,len(b.triangles),3):
#    print i/3, b.triangles[i]-1,b.triangles[i+1]-1,b.triangles[i+2]-1
print area(b)/100.
#print b.numberoftriangles
print Ix(b)/10000.
# plot_triangles(b)
#trigen.vtk("test.vtk", b)

